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We determine the phase diagram of liquid supercooled water by extensive computer simulations 
using the TIP5P-E model [J. Chem. Phys. 120, 6085 (2004)]. We find that the transformation of 
water into a low density liquid in the supercooled range strongly enhances the solubility of hydropho- 
bic particles. The transformation of water into a tetrahedrally structured liquid is accompanied by a 
minimum in the hydration entropy and enthalpy. The corresponding change in sign of the solvation 
heat capacity indicates a loss of one characteristic signature of hydrophobic hydration. The observed 
behavior is found to be qualitatively in accordance with the predictions of the information theory 
model of Garde et al. [Phys. Rev. Lett. 77, 4966 (1996)]. 



The thermodynamical anomalies of liquid water are 
considered to be caused by a transformation between two 
different liquid forms of water buried in the deeply super- 
cooled region The two differently dense liquids have 
well characterized counterparts in the glassy state: The 
(very) high density and low density amorphous ice forms 
0, Eg • Computer simulation studies have furnished a pic- 
ture of a first order liquid-liquid phase transition between 
two liquids ending up in a metastable critical point 0,0- 
Although singularity free scenarios might as well explain 
the properties of supercooled water |l(, there is experi- 
mental support for the liquid-liquid critical point hypoth- 
esis from the changing slope of the metastable melting 
curves observed for different ice polymorphs 0, • To 
make the situation even more puzzling, recent computer 
simulations provide evidence that there might be even 
more than one liquid- liquid transition ||. 

One prominent anomaly of liquid water is the increas- 
ing solubility of hydrophobic gases with decreasing tem- 
perature . This behavior is the consequence of a nega- 
tive solvation entropy of small hydrophobic particles [lOj . 
Since the corresponding solvation enthalpy is also nega- 
tive, the observed low solubility (the large positive sol- 
vation free energy) of small hydrophobic particles is due 
to a dominating entropy effect. The origin of the nega- 
tive hydration entropy is widely regarded as being due 
to the bias in the hydration waters orientational space, 
as the water molecules are trying to preserve their hy- 
drogen bond network [ill ]. Entropy and enthalpy effects 
have been shown to be determined by the water in the 
first hydration shell I n addition, the broken hydro- 
gen bond states, with hydrogen bond donors/acceptors 
pointing towards the hydrophobic particles, are increas- 
ingly populated with increasing temperature 0. Hence, 
the positive solvation heat capacity of hydrophobic par- 
ticles is considered as one of the key signatures for hy- 
drophobic hydration 0]. A recent experimental study 
of Souda investigating alkane layers adsorbed on an 
amorphous solid water substrate, shows that the alkane 
phase gets soaked into the water phase in the tempera- 
ture region close to the suspected glass transition tem- 



perature (T^j. One possible explanation (among others) 
might be that the highly viscous low density liquid form 
of water provides a significantly increased solubility for 
hydrophobic molecules. The scope of the present contri- 
bution is to monitor hydration of a small hydrophobic 
particle as we penetrate the deeply supercooled region. 
Our particular interest is to elucidate how the transfor- 
mation into a highly tetrahedrally ordered low density 
liquid (LDL) phase [lfj affects hydrophobic hydration. 

An important problem regarding simulations of highly 
viscous liquids close to the glass transition is to pro- 
vide proper sampling. To overcome this problem we 
perform parallel tempering simulations of an extended 
ensemble of states |l7j using the technique of volume- 
temperature replica exchange molecular dynamics simu- 
lation (VTREMD) [13. To represent liquid supercooled 
water, we employ the recently proposed TIP5P-E model 
for water, which was fitted to reproduce waters density 
maximum and related thermodynamical anomalies while 
treating the Coulomb interactions with the Ewald sum 
[raj . For the VTREMD simulation [H we consider a 
grid of 440 (V,T)-states 21]. Starting from a set of equi- 
librated initial configurations obtained at ambient con- 
ditions, the VTREMD simulation |22J was conducted for 
20 ns, providing a total 8.8 /xs worth of trajectory data. 
The average time interval between two successful state- 
exchanges was obtained to be about 3ps. During the 
entire course of the simulation each replica has crossed 
the whole temperature and density interval several times. 
After an initial equilibration period of about 4 ns, the av- 
erage pressure and potential energies show convergence 
even for the lowest temperatures and the remaining 16 ns 
are used for analysis. 

Figure 1 shows the phase diagram of liquid supercooled 
water in terms of a contour plot of the P(p, T)-data as 
obtained by the VTREMD simulations. The TIP5P-E 
phase diagram apparently exhibits a first order phase 
transition between two metastable liquid phases, ending 
in a second critical point C* . The location of the coexis- 
tence line, estimated from a Maxwell construction using 
the subcritical isotherms shown in Figure |2Ji, has to be 
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FIG. 1: Phase diagram for stable and supercooled liquid 
TIP5P-E water. The spacing between different contour 
colours corresponds to a pressure drop of 25 MPa. Selected 
isobars are indicated. The heavy dashed line denotes the high 
density / low density liquid (HDL/LDL) coexistence line. The 
HDL/LDL critical point is located at C* with T* = 210 K, 
P* =310 MPa and p* = 1.09gcm -3 . Darker shading indicates 
the LDL-basin. 

seen as a rough guess only, since the exact form of the van 
der Waals loops in the two phase region might depend on 
the system size. The location of C* is close the values 
reported by Yamada et al. for the original TIP5P model 
with T* = 217K, P*=340MPaand p* = 1.13 gem" 3 0. 
Figure 2a compares experimental data [24| with several 
selected isobars according to a linear interpolation from 
the TIP5P-E P(p, T) data-set. Although the location of 
the density maximum at normal pressure is close to the 
experimental values, the TIP5P-E model exhibits a con- 
siderably larger thermal expansivity as real water and is 
significantly more compressible. Figure 2b shows that the 
isobars derived from VTREMD data-set match exactly 
the data obtained from conventional A^PT-simulations. 
For comparison the isobar according to Rick's simulation 
is shown (lij. We denote small, but significant differ- 
ences. Since both simulations were using Ewald summa- 
tion, their origin is not quite clear at present. We would 
like to emphasize one particular detail of the low pressure 
isobars, indicated in Figures 1 and 2b: The density goes 
through a minimum after the transformation into the low 
density liquid has taken place. Apparently at low tem- 
peratures the expansivity seems to behave as it would be 
expected for a conventional liq uid. I. Brovchenko et al. 
H and recently P. H. Poole |25j have made similar obser- 
vations for the ST2-model. Figure 2c shows the evolution 
of the composite radial distribution function h(r) along 
the 222.7 K-isotherm. The highlighted patterns obtained 
for the lowest and highest densities show strong simi- 
larity with the functions determined experimentally by 
Bellissent-Funel et al. for high and low density amor- 




FIG. 2: (a): Selected isobars for TIP5P-E water (dot- 
ted lines) compared with experimental data 24]. (b) 
Two isobars for TIP5P-E. Open circles: Interpolated 
data from the VTREMD P(p,T) data-set. Open 
squares: Constant temperature/pressure (NPT) MD 
simulations. Filled circles: Ref. fl^|. (c) h(r) = 
Anr{N/V) [0.092 g OQ (r) +0.486 ff HH (r) + 0.422 g QR (r) - l]- 
functions for the 222.7 K isotherm. The functions are shifted 
by an increment of 10nm~ 2 . (d) Near- and sub-critical 
isotherms. 

phous ice 26] . The observed low density pattern has 
been shown to be related to a high tetrahedral order of 
water's first hydration shell ^(|. Similarities between the 
local structure of LDA and ice Ih have also been demon- 
strated recently by Finney et al. from the analysis of 
neutron scattering experiments [22} • We do not observe 
indications for the presence of further liquid-liquid tran- 
sitions for TIP5P-E, as suggested by Brovchenko et al. 
for the ST2 model [8j . However, further phase transitions 
might be present at even lower temperatures. Finally, no 
signs for crystallization are observed, as denoted by Ya- 
mada et al. 0. Probably our smaller system size and 
the random walk of the replicas through state-space, reg- 
ularly reaching stable regions, prevents the systems from 
crystallizing. 

The hydrophobic hydration of a Lennard- Jones Ar gon 
particle (u^^q = 0.3290 nm and e^ Q/fc = 98.9 K) |30| 
is given as excess chemical potential p ex — — In 7 for 
infinite dilution, where 7 is the solubility and 13 = 1/ kT. 
We employ the Widom particle insertion method 0] 
with p ex = — (exp(— /?$)), where $ is the energy 
of an inserted test-particle. The brackets (...) indicate 
canonical sampling. Details about the calculation are 
given elsewhere 30]. Figure 3 shows the excess chemical 
potential for Argon dissolved in TIP5P-E water for all 
simulated state points. First of all we would like to em- 
phasize that the solubility increases strongly when pene- 
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FIG. 3: Excess chemical potential of Ar dissolved in TIP5P- 
E water as a function of P and T. Thin lines indicate all 
calculated isotherms ranging from 170 K to 356.3 K. The in- 
sert shows experimental data for 298 K |2fl| and simulated 
isotherms close-by. Thick lines indicate the LDL-basin and 
the HDL/LDL-coexistence region (see Figure Ql. 
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FIG. 4: The symbols indicate six isobars, interpolated from 
the /iea; (p, T) TIP5P-E data-set. The dashed line indicates ex- 
perimental data for 0.1 MPa |29|. Full lines represent predic- 
tions for p e;E = ap 2 T+b with a = 8xl0~ 2 kj mol -1 K -1 cm 6 g~ 2 
and b = — 16 kj mol -1 . The insert shows enthalpy (h ex ) and 
entropy (—Ts ex ) contributions to fi ex for the OMPa isobar. 



trating the deeply supercooled region with p ex becoming 
even negative. The insert of Figure|3]compares the exper- 
imental data of Kennan and Pollack [2^] with data from 
our simulation. The pressure dependence is expressed by 
a positive partial molar volume of 30 cm 3 mol -1 , close 
to the experimental value of 25 cm mol" 1 [HI. The 
decreasing solubility upon density increase basically re- 
flects the loss of available free volume. A larger par- 
tial molar volume in the simulation is hence in line with 
the increased isothermal compressibility observed for the 
TIP5P-E model, indicated by Figure^. Two regimes of 
small and strong pressure dependence are denoted, obvi- 
ously related to the transformation between the high and 



low density liquid forms of water. We would like to point 
out that we find an increase of the partial molar vol- 
ume for Ar around normal pressure for the 260 — 255 K 
isotherms. Kennan and Pollack made a similar obser- 
vation for Xenon around 298 K. If, and how these two 
observations are related has to be further investigated. 
Figure 4 shows the temperature dependence of p ex for 
several selected isobars ranging from —200 MPa to the 
supercritical 350 MPa. The experimental 0.1 MPa isobar 
[2jj is given for comparison. From computer simulations 
of water, S. Garde et al. have derived an information 
theory (IT) model [13, H3, giving simple analytic ex- 
pressions for the hydrophobic hydration as a function of 
temperature and density. The leading term in the IT 
model strongly suggests a quadratic relation between the 
excess chemical potential and the water number density 
according to p ex /k p 2 Tv 2 /2a 2 [33|, where v denotes 
the volume of a hydrophobic hard sphere particle, while 
ofj = (n 2 ) — (n) 2 indicates the variance of the number of 
water molecules in a sphere of volume v. The lines in Fig- 
ure 4 indicate a temperature dependence as suggested by 
the IT model, assuming the term v 2 /2a^ to be constant 
and shifting the isobars by a constant offset to account 
for attractive interactions. Well reproduced is the change 
in slope when passing the transformation into the low 
density liquid, qualitatively correct even for the super- 
critical isotherm at 350 MPa, suggesting a discontinuity 
of fi ex . Hence an important conclusion is that the weak 
monotonous temperature and density dependence of fi ex , 
as proposed by the IT model, and confirmed for the high 
temperature regime l32i . is apparently also valid for the 
deeply supercooled region. This seems to be in line with 
the continuous structural transformation of water as a 
function of density (Figure [2p) and temperature (not 
shown) . The thermodynamics of hydrophobic hydration 
of small solutes is therefore largely predetermined by the 
location of the water isobars in the p, T-plane. As a con- 
sequence, also a waterlike liquid showing a continuous 
singularity free transition might exhibit a quite similar 
//^-dependence (apart from the discontinuous behavior 
of fi ex above P c ). As shown in the insert of Figure 4, the 
change in slope of p ex is related to extrema in the temper- 
ature dependence of the entropy and enthalpy contribu- 
tions. The minimum in the enthalpy causes a change in 
sign of the corresponding heat capacity from positive at 
high temperatures to negative for the low density liquid. 
Apparently, a prominent signature of the "hydrophobic 
hydration" vanishes while passing the transition at about 
250 K. In addition, an increasing slope of h ex upon ap- 
proaching 250 K, as seen in the insert, and more promi- 
nently pronounced with increasing pressure (not shown) , 
is consistent with the experimentally observed increase of 
the solvation heat capacity upon cooling [13J ■ Following 
the interpretation of the hydrophobic hydration accord- 
ing to Dill et al. 11], the change in sign of C p is explained 
as follows: In the LDL-regime, water's coordination num- 
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ber approaches four, while the structure corresponds to a 
continuous random tetrahedral network. Due to the lack 
of additional possible hydrogen binding partners, the en- 
vironment of a water molecule in the LDL-bulk and in the 
hydrophobic hydration shell has greater similarity than 
in the "high density liquid" at ambient temperatures. 
This structural convergence in bulk and hydration shell 
in the LD-liquid leads to a smaller entropy difference and 
hence to an increase in s ex . Moreover, a larger hydropho- 
bic particle in the LD-liquid might even introduce disor- 
der and hydrogen bond breaking into the (over) stretched 
hydrogen bond network and thus might even experience 
a positive hydration entropy. A further investigation on 
particle size dependence, hydrophobic aggregates and hy- 
drophobic interaction might reveal a behavior markedly 
different from that at ambient conditions. 
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